{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Dynamic mode decomposition with control on a 2D linear system \n",
    "\n",
    "Dynamic mode decomposition with control (DMDc) aims to disambiguate the effect of control/actuation\n",
    "from the unforced dynamics.\n",
    "We apply DMDc to a low-dimensional, linear system\n",
    "(this is example in Sec. 4 in Proctor et al., _\"Dynamic Mode Decomposition with Control\"_, SIAM 2016).\n",
    "\n",
    "$$x_{k+1} =\\begin{bmatrix} 1.5 & 0\\\\ 0 & 0.1 \\end{bmatrix}x_k + \\begin{bmatrix} 1\\\\ 0 \\end{bmatrix} u_k$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We first import the pyKoopman package and other packages for plotting and matrix manipulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../src')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "import pykoopman as pk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define state and control matrices of the linear control system and collect data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X1 =  [[4.000000e+00 7.000000e+00]\n",
      " [2.000000e+00 7.000000e-01]\n",
      " [1.000000e+00 7.000000e-02]\n",
      " [5.000000e-01 7.000000e-03]\n",
      " [2.500000e-01 7.000000e-04]\n",
      " [3.750000e-01 7.000000e-05]\n",
      " [1.062500e+00 7.000000e-06]\n",
      " [2.593750e+00 7.000000e-07]\n",
      " [6.890625e+00 7.000000e-08]]\n",
      "X2 =  [[2.00000000e+00 7.00000000e-01]\n",
      " [1.00000000e+00 7.00000000e-02]\n",
      " [5.00000000e-01 7.00000000e-03]\n",
      " [2.50000000e-01 7.00000000e-04]\n",
      " [3.75000000e-01 7.00000000e-05]\n",
      " [1.06250000e+00 7.00000000e-06]\n",
      " [2.59375000e+00 7.00000000e-07]\n",
      " [6.89062500e+00 7.00000000e-08]\n",
      " [1.53359375e+01 7.00000000e-09]]\n",
      "C =  [[-4. ]\n",
      " [-2. ]\n",
      " [-1. ]\n",
      " [-0.5]\n",
      " [ 0. ]\n",
      " [ 0.5]\n",
      " [ 1. ]\n",
      " [ 3. ]\n",
      " [ 5. ]]\n"
     ]
    }
   ],
   "source": [
    "from pykoopman.common  import advance_linear_system\n",
    "\n",
    "A = np.array([[1.5, 0],[0, 0.1]])\n",
    "B = np.array([[1],[0]])\n",
    "\n",
    "x0 = np.array([4,7])\n",
    "u = np.array([-4, -2, -1, -0.5, 0, 0.5, 1, 3, 5])\n",
    "n = len(u)+1\n",
    "x,_ = advance_linear_system(x0,u,n,A,B)\n",
    "X1 = x[:-1,:]\n",
    "X2 = x[1:,:]\n",
    "C = u[:,np.newaxis]\n",
    "print('X1 = ', X1)\n",
    "print('X2 = ', X2)\n",
    "print('C = ', C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Apply first the standard DMD to the state data collected from the controlled system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[ 2.17160989e+00 -9.95420579e-01]\n",
      " [-1.58023594e-17  1.00000000e-01]]\n"
     ]
    }
   ],
   "source": [
    "U, s, Vh = np.linalg.svd(X1.T, full_matrices=False)\n",
    "Aest = np.dot(X2.T,np.dot(Vh.T*(s**(-1)),U.T))\n",
    "print('A = ', Aest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is obviously not correct.\n",
    "So let's apply DMDc on the data from the controlled system.\n",
    "We assume for now that the control matrix B is known."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>Koopman(observables=Identity(),\n",
       "        regressor=DMDc(input_control_matrix=array([[1],\n",
       "       [0]]),\n",
       "                       svd_output_rank=3, svd_rank=3))</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" ><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">Koopman</label><div class=\"sk-toggleable__content\"><pre>Koopman(observables=Identity(),\n",
       "        regressor=DMDc(input_control_matrix=array([[1],\n",
       "       [0]]),\n",
       "                       svd_output_rank=3, svd_rank=3))</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" ><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">observables: Identity</label><div class=\"sk-toggleable__content\"><pre>Identity()</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-3\" type=\"checkbox\" ><label for=\"sk-estimator-id-3\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">Identity</label><div class=\"sk-toggleable__content\"><pre>Identity()</pre></div></div></div></div></div></div><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-4\" type=\"checkbox\" ><label for=\"sk-estimator-id-4\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">regressor: DMDc</label><div class=\"sk-toggleable__content\"><pre>DMDc(input_control_matrix=array([[1],\n",
       "       [0]]), svd_output_rank=3,\n",
       "     svd_rank=3)</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-5\" type=\"checkbox\" ><label for=\"sk-estimator-id-5\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">DMDc</label><div class=\"sk-toggleable__content\"><pre>DMDc(input_control_matrix=array([[1],\n",
       "       [0]]), svd_output_rank=3,\n",
       "     svd_rank=3)</pre></div></div></div></div></div></div></div></div></div></div>"
      ],
      "text/plain": [
       "Koopman(observables=Identity(),\n",
       "        regressor=DMDc(input_control_matrix=array([[1],\n",
       "       [0]]),\n",
       "                       svd_output_rank=3, svd_rank=3))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DMDc = pk.regression.DMDc(svd_rank=3, input_control_matrix=B)\n",
    "model = pk.Koopman(regressor=DMDc)\n",
    "model.fit(x,u=C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.09412963 0.63520687]\n",
      " [0.63520687 0.50587037]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aest = model.A[:,:2]\n",
    "Best = model.A[:,-1:]\n",
    "\n",
    "print(Aest)\n",
    "np.allclose(A,Aest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "This yields the correct system matrix. Let's further assume B is also unknown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "DMDc = pk.regression.DMDc(svd_rank=3)\n",
    "\n",
    "model = pk.Koopman(regressor=DMDc)\n",
    "model.fit(x,u=C)\n",
    "Aest = model.ur @ model.A @ model.ur.T\n",
    "\n",
    "Best = model.ur @ model.B\n",
    "\n",
    "assert np.allclose(B,Best)\n",
    "assert np.allclose(A,Aest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we can simulate the system using the learned DMDc model\n",
    "and compare with the true solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x22f8cc7ed40>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/IAAAEqCAYAAAC2kGtTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiU0lEQVR4nO3deVxU9f7H8dcwyKrggoooivsu7gpaWlpmZtmiZZaaZrfS0rxtdsu0Rdqzm161cmsxK29Zvxa7aqWZ+0LuiguiKSqmIKggzPn9MTKKAoLOcGaG9/PxOA9mzpwz854R+c7nnO/5fi2GYRiIiIiIiIiIiEfwMTuAiIiIiIiIiBSdCnkRERERERERD6JCXkRERERERMSDqJAXERERERER8SAq5EVEREREREQ8iAp5EREREREREQ+iQl5ERERERETEg6iQFxEREREREfEgvmYHcEc2m42DBw9Srlw5LBaL2XFEREQwDIOTJ08SERGBj4+OwzuD2nsREXEnxWnrVcjn4+DBg0RGRpodQ0RE5BL79++nRo0aZsfwCmrvRUTEHRWlrVchn49y5coB9g8wJCTE5DQiIiKQlpZGZGSko42Sq6f2XkRE3Elx2noV8vnI7V4XEhKihl1ERC4rKQlSUoq/X1gY1KxZvH3UBdx5rqS9L8l/axERKZ2K0tarkBcREbkKSUnQsCGcOVP8fQMCYMcOFXieQv/WIiLiLjRajoiIyFVISbmywg7s+13J2V0xh/6tRUTEXaiQFxEREREREfEg6lp/hQzDIDs7m5ycHLOjeK0yZcpgtVrNjiEickXqkcAQZhBFIolEMYMh7KK+2bG8RlRUFPv27btk/aOPPsrkyZNNSCQiIlJyVMhfgaysLA4dOsSpU6fMjuLVLBYLNWrUoGzZsmZHEREplsHM5CMexMCCBQMDC0/zBkOZzmwGmx3PK6xZsybPwfTNmzdzww030LdvXxNTiYiIlAwV8sVks9nYu3cvVquViIgI/Pz8NIKwCxiGwdGjRzlw4AD169fXmXkR8Rj1SOAjHsSKLc96A5jOUJbRmd3UMyecF6lcuXKe+6+99hp169alS5cuJiUSEREpOSrkiykrKwubzUZkZCRBQUFmx/FqlStXJjExkbNnz6qQFxGPMYQZGFx6gNcCGFgYynSeI67kg3mxrKwsPv30U0aPHl3owfXMzEwyMzMd99PS0koinoiIiNOpkL9CPj7FGydQ884Wn3o6iIgniiIRC0YBjxpEkViScUqF+fPnc+LECQYPHlzodnFxcYwfP75kQomIiLiQ6aPWL126lN69exMREYHFYmH+/Pl5Hh88eDAWiyXPctNNN132eSdPnkxUVBQBAQF06NCB1atXu+gdXF7uvLNt2hR/adjQvr+IiHiGRKLyPSNvZyGRqJKMUypMnz6dnj17EhERUeh2Y8aMITU11bHs37+/hBKKiIg4l+mFfEZGBtHR0YWOMHvTTTdx6NAhx/L5558X+pxffPEFo0eP5sUXX2T9+vVER0fTo0cPjhw54uz4RaJ5Z0VESo8ZDDk3wF1eBmDBYDpDzYjltfbt28eiRYt48MEHL7utv78/ISEheRYRERFPZHoh37NnT1555RVuv/32Arfx9/cnPDzcsVSoUKHQ53znnXcYNmwYDzzwAE2aNGHq1KkEBQUxY8YMZ8cXERHJYxf1+YBhAOTgQzY+ZGPFhg9Dma6B7pxs5syZVKlShV69epkdRUREpMSYXsgXxW+//UaVKlVo2LAhjzzyCMeOHStw26ysLNatW0f37t0d63x8fOjevTsrVqzId5/MzEzS0tLyLN7m4ssTLl7GjRtndkQREa9Rlz1YgOXE8hX9eJOnaMgOTT3nZDabjZkzZzJo0CB8fTXsj4iIlB5u3+rddNNN3HHHHdSuXZvdu3fz3HPP0bNnT1asWJHvSOYpKSnk5ORQtWrVPOurVq3K9u3b832N0jD4zaFDhxy3v/jiC8aOHcuOHTsc6y6cq90wDHJycvSlSETkClTjIN1YDMBgZrGHuiYn8l6LFi0iKSmJIUOGmB1FRESkRLn9Gfl77rmHW2+9lebNm9OnTx++//571qxZw2+//ea017jawW8MAzIyCl5On766fKdPF/zcRkEDI1/kwksTQkNDsVgsjvvbt2+nXLly/PTTT7Rp0wZ/f3+WLVvG4MGD6dOnT57nGTVqFF27dnXct9lsxMXFUbt2bQIDA4mOjmbevHlX94ZFRDzYvczBio0/iFUR72I33ngjhmHQoEEDs6OIiIiUKI875VqnTh3CwsLYtWsX3bp1u+TxsLAwrFYrhw8fzrP+8OHDhIeH5/uc/v7++Pv7X3GmU6fgghPaTte5c8GPpadDcLBzXufZZ5/lrbfeok6dOpcdhyBXXFwcn376KVOnTqV+/fosXbqU++67j8qVK9OlSxfnBBMRcWNhYRAQcH5Q0/v5BIBPuP+y+wYE2PcXz3Dxv3Vx6N9aREScyeMK+QMHDnDs2DGqVauW7+N+fn60adOGxYsXO84m22w2Fi9ezIgRI0owqed56aWXuOGGG4q8fWZmJhMmTGDRokXExMQA9gMty5YtY9q0aSrkRaRUqFkTduywzzByZNFGop/ZSCZ+3PdtPx6qUfi+YWH2/cUzXPhvnZ8//oDHH4fISLhoNl39W4uIiFOZXsinp6eza9cux/29e/cSHx9PxYoVqVixIuPHj+fOO+8kPDyc3bt38/TTT1OvXj169Ojh2Kdbt27cfvvtjkJ99OjRDBo0iLZt29K+fXsmTpxIRkYGDzzwgEveQ1CQ/cx4QeLjCz+rfjnLlkHLlgW/trO0bdu2WNvv2rWLU6dOXVL8Z2Vl0apVK+cFExFxczVr2pffHv4UgA0Rveh8a0WTU4kr5P5b56d2bXshv3+/vZivXLlks4mISOlheiG/du1arrvuOsf90aNHAzBo0CCmTJnCxo0bmT17NidOnCAiIoIbb7yRl19+OU9X+N27d5NyweHxu+++m6NHjzJ27FiSk5Np2bIlCxYsuGQAPGexWArv3h4YeHXPHxjovO7zhQm+6EV8fHwwLroI/+zZs47b6eeOXvzwww9Ur149z3ZXc6mCiIgnysmBrxJa4kcM1oEDzY4jJqhQARo3hm3bYMUKuPVWsxOJiIi3Mr2Q79q16yXF4oV+/vnnyz5HYmLiJetGjBihrvRXqXLlymzevDnPuvj4eMqUKQNAkyZN8Pf3JykpSd3oRaTU++UX+M+Je5lb8V4OjSviSKTidWJj7YX88uUq5EVExHXcftR6Mc/111/P2rVr+fjjj0lISODFF1/MU9iXK1eOJ598kieeeILZs2eze/du1q9fz/vvv8/s2bNNTC4iUvI+sY9xx913g5+/xdwwYprYWPvP5cvNzSEiIt5NhbwUqEePHrzwwgs8/fTTtGvXjpMnTzLwou6iL7/8Mi+88AJxcXE0btyYm266iR9++IHatWublFpEpORlHMmgwhdTqcDf3H/5werFi+UW8mvWQFaWuVlERMR7WYzC+rWXUmlpaYSGhpKamkpISEiex86cOcPevXupXbs2AQEBRXq+9euhTZsrz7NuHbRufeX7e6or+axFRMzwxyOf0mnq/ewq05i6mVuxuOCEfGFtk1wZV3ymNpt9kLu//4bVq6FdO6c8rYiIlALFaZd0Rr4E5M47eyU076yIiPvz/8rer/5AbD+XFPHiOXx8oGNH+211rxcREVcxfbC70uBy884WRvPOioi4t8Pxh2h1bBEAtV+4z+Q04g5iY+HHH+2F/MiRZqcRERFvpEK+hBQ276yIiHiubWM/pys2NpWLoXm3embHETegAe9ERMTV1LVeRETkKoQvtHerP95Lo9yJXbt2YLXCgQOwf7/ZaURExBupkBcREblCO7/eTKMz8WRRhuYv9zM7jriJsmUhOtp+W2flRUTEFVTIi4iIXKH4D1eTjZX11XpRoV4ls+OIG8ntXr9ihbk5RETEO6mQFxERuQI5OfDExiFU5y9Sn3vD7DjiZmJi7D91Rl5ERFxBhbyIiMgV+PVXOHgQzlaoStdh9c2OI24m94z8hg1w6pS5WURExPuokBcREbkC/51+AoB+/cDf39ws4n5q1YJq1SA7G9auNTuNiIh4GxXypcjgwYOxWCxYLBbKlClD1apVueGGG5gxYwY2m82xXVRUFBaLhblz517yHE2bNsVisTBr1qxLtrdYLAQGBhIVFUW/fv345ZdfSuJtiYiUuIwjGbw+tya/0pUHeqeYHUfckMWiaehERMR1VMibKSEBxoyB/v3tPxMSXP6SN910E4cOHSIxMZGffvqJ6667jpEjR3LLLbeQnZ3t2C4yMpKZM2fm2XflypUkJycTHBx8yfO+9NJLHDp0iB07dvDxxx9Tvnx5unfvzquvvury9yQiUtLix39LCCep7buf9j01yJ3kT4W8iIi4igp5s8ycCY0awZtvwpdf2n82agQXnOl2BX9/f8LDw6levTqtW7fmueee49tvv+Wnn37Kc5Z9wIABLFmyhP0XTIA7Y8YMBgwYgK+v7yXPW65cOcLDw6lZsybXXnstH3zwAS+88AJjx45lx44dju22bNnCLbfcQkhICOXKleOaa65h9+7dLn3PIiLO5velfe74vbH3YfGxmJym9Prrr7+47777qFSpEoGBgTRv3py1btSP/cKR6w3D3CwiIuJdVMg7U0ZGwcuZM+e3S0iABx8Em80+7PGFP4cOhc2bi/a8TnL99dcTHR3N119/7VhXtWpVevTowezZswE4deoUX3zxBUOGDCny844cORLDMPj2228B+xeua6+9Fn9/f3755RfWrVvHkCFD8vQEEBFxd4f/TKZ1yv8AiHr+PpPTlF7Hjx+nU6dOlClThp9++omtW7fy9ttvU6FCBbOjObRqBX5+kJICu3aZnUZERLzJpadW5cqVLVvwYzffDD/8YL89Y4a9aM+PzWbfNinp/LqoKPu3gIs58fB+o0aN2LhxY551Q4YM4Z///Cf/+te/mDdvHnXr1qVly5ZFfs6KFStSpUoVEhMTAZg8eTKhoaHMnTuXMmXKANCgQQNnvQURkRKxbezndMXGprIdaX6DRqs3y+uvv37JZWC1a9c2MdGl/P2hbVt71/rly6G+fl1ERMRJdEbeDOcK2wJdePa+hBiGgcWSt3tor169SE9PZ+nSpcyYMaNYZ+Pze974+HiuueYaRxEvIuKJqi60d6v/u9f9Jicp3b777jvatm1L3759qVKlCq1ateLDDz8sdJ/MzEzS0tLyLK6m6+RFRMQVdEbemdLTC37Maj1/OyrKfj8nJ//tBg3Ku+5yhb8TbNu27ZIzGb6+vtx///28+OKLrFq1im+++aZYz3ns2DGOHj3qeN7AwECn5RURMUPC/C00Pr2BLMrQ/OW7zY5Tqu3Zs4cpU6YwevRonnvuOdasWcPjjz+On58fgy5uR8+Ji4tj/PjxJZpThbyIiLiCzsg7U3BwwUtAwPnthgwpuFu8YcA//lG053WSX375hU2bNnHnnXde8tiQIUNYsmQJt912W7GvO3zvvffw8fGhT58+ALRo0YLff/+ds2fPOiO2iEiJm70kivv5mHmNnqdifY1WbyabzUbr1q2ZMGECrVq14qGHHmLYsGFMnTq1wH3GjBlDamqqY7lwQFdXiYmx/9yyBVJTXf5yIiJSSqiQN0P9+jB9Ovj42M/AX/hz+nSoV89lL52ZmUlycjJ//fUX69evZ8KECdx2223ccsstDBw48JLtGzduTEpKyiVT0V3s5MmTJCcns3//fpYuXcpDDz3EK6+8wquvvkq9c+9nxIgRpKWlcc8997B27VoSEhL45JNP8oxqLyLirnJyYNZXwXzK/fi/MtbsOKVetWrVaNKkSZ51jRs3JunCMWYu4u/vT0hISJ7F1cLDoU4d+3H6Vatc/nIiIlJKqGu9WQYPhs6d7YV7YqK9u/3QoS4t4gEWLFhAtWrV8PX1pUKFCkRHR/Pvf/+bQYMG4eOT/3GdSpUuf9Zp7NixjB07Fj8/P8LDw+nYsSOLFy/muuuuy/M8v/zyC0899RRdunTBarXSsmVLOnXq5LT3JyLiKr/9Bn/9BeXLwy23mJ1GOnXqdMmB4J07d1KrVi2TEhUsNhb27LF3r7/xRrPTiIiIN1Ahb6Z69SAursRebtasWXnmii9I4mWuyT9x4kSxtr9QixYt+Pnnn4u8vYiIu0h96mWeJACj9yD8/auYHafUe+KJJ4iNjWXChAn069eP1atX88EHH/DBBx+YHe0SMTHw6ae6Tl5ERJzH9K71S5cupXfv3kRERGCxWJg/f77jsbNnz/LMM8/QvHlzgoODiYiIYODAgRw8eLDQ5xw3bhwWiyXP0qhRIxe/ExER8VanUk5xw4Y3eJOnGdhBlwO5g3bt2vHNN9/w+eef06xZM15++WUmTpzIgAEDzI52idwB71auzH+cWxERkeIy/Yx8RkYG0dHRDBkyhDvuuCPPY6dOnWL9+vW88MILREdHc/z4cUaOHMmtt97K2rVrC33epk2bsmjRIsd9X1/T36qIiHio+PHfEks6+32jaP6wLgdyF7fccgu3eMB1Ds2aQdmycPKkfdC7Fi3MTiQiIp7O9Oq2Z8+e9OzZM9/HQkNDWbhwYZ51kyZNon379iQlJVGzZs0Cn9fX15fw8HCnZhURkdKpzBf2ueN3d7yPSKvpndnEw/j6QocOsHixvXu9CnkREblaHvdtJDU1FYvFQvny5QvdLiEhgYiICOrUqcOAAQMKHcVWRESkIEc3H6bV0f8BEPXC/SanEU+V271+xQpzc4iIiHfwqEL+zJkzPPPMM/Tv37/QKWM6dOjArFmzWLBgAVOmTGHv3r1cc801nDx5Mt/tMzMzSUtLy7OIiIgAbH3hc3zJYUtwe6JubGB2HPFQuYW8BrwTERFn8JhC/uzZs/Tr1w/DMJgyZUqh2/bs2ZO+ffvSokULevTowY8//siJEyf48ssv890+Li6O0NBQxxIZGXnZPIZhXNH7kKLTZywi7qDKz/Zu9cdu1tl4uXIdOth/7toFR46Ym0VERDyfRxTyuUX8vn37WLhwYaFn4/NTvnx5GjRowK5du/J9fMyYMaSmpjqW/fv3F/hcZcqUAewD8YlrZWVlAWC1Wk1OIiKl1bYNZ/jzdH3SKEfTl+8xO454sAoVoEkT+211rxcRkatl+mB3l5NbxCckJPDrr79SqVKlYj9Heno6u3fv5v778z+b4u/vj7+/f5Gey2q1Ur58eY6cO5weFBSExWIpdiYpnM1m4+jRowQFBWnGARExzcdfBvAac7mz1xnmNQwwO454uNhY2LrV3r3+ttvMTiMiIp7M9AopPT09z5nyvXv3Eh8fT8WKFalWrRp33XUX69ev5/vvvycnJ4fk5GQAKlasiJ+fHwDdunXj9ttvZ8SIEQA8+eST9O7dm1q1anHw4EFefPFFrFYr/fv3d0rm3NHwj6hvnEv5+PhQs2ZNHSgREVPYbPDZZ/bb9wxWES9XLzYWPvpI18mLiMjVM72QX7t2Ldddd53j/ujRowEYNGgQ48aN47vvvgOgZcuWefb79ddf6dq1KwC7d+8mJSXF8diBAwfo378/x44do3LlynTu3JmVK1dSuXJlp2S2WCxUq1aNKlWqcPbsWac8p1zKz88PHx+PuPpDRLzQ6s8SCN6fTWhoYzxgqnLxALkD3q1dC1lZcO58hIiISLGZXsh37dq10EHNijLgWWJiYp77c+fOvdpYRWK1WnX9toiIlzr7chzbmMm39cYTEDDW7DjiBRo0gIoV4e+/IT4e2rc3O5GIiHgqne4UERG5yKmUU0QnzAMgalBXc8OI17BYICbGflvd60VE5GqokBcREblI/EvfEcJJDlhr0fyRzmbHES+i+eRFRMQZVMiLiIhcpMxc+9zxuzreh4+vmkpxntxC/o8/oAhXD4qIiORL305EREQucHTzYVod/RmAWs/nP22pyJVq1w6sVjh4EPbvNzuNiIh4KhXyIiIiF9g6di6+5LAluB21b2podhzxMsHBkDsRz4oVpkYREREPpkJeRETkAgGL/g+AlJ4DTU4i3krXyYuIyNVSIS8iInLOtm1w7ckfuNPnG5q9co/ZccRLaeR6ERG5WirkRUREzvnkE8jCn7O9+lCpYZjZccRL5Z6R37ABMjLMzSIiIp5JhbyIiAhgyzH47FP7MOL3a4w7caGaNSEiAnJyYO1as9OIiIgnUiEvIiICbHx/Cb/ur8vzAW/Su7fZacSbWSy6Tl5ERK6OCnkRERHg5H8+oQ57ubHmDgICzE4jlzNu3DgsFkuepVGjRmbHKrLcQl4j14uIyJXwNTuAiIiI2U7/fZoWCfMACBmufvWeomnTpixatMhx39fXc77WXHhG3jDsZ+lFRESKynNaPBERERfZMP47YknjgLUmzR+9xuw4UkS+vr6Eh4ebHeOKtGoF/v5w7BgkJECDBmYnEhERT6Ku9SIiUur5zv0EgF0d7sPHV02jp0hISCAiIoI6deowYMAAkpKSCt0+MzOTtLS0PItZ/PygbVv7bV0nLyIixaVvKyIiUqod3XKE1kcWAFDreXWr9xQdOnRg1qxZLFiwgClTprB3716uueYaTp48WeA+cXFxhIaGOpbIyMgSTHwpDXgnIiJXSoW8iIiUalte/AJfctga1JbaPT1nsLTSrmfPnvTt25cWLVrQo0cPfvzxR06cOMGXX35Z4D5jxowhNTXVsezfv78EE19KhbyIiFwpXSMvIiKl2vStsexiKLVvi6GJ2WHkipUvX54GDRqwa9euArfx9/fH39+/BFMVLibG/nPrVjhxAsqXNzONiIh4Ep2RFxGRUmv7dvh0Wxse8f2IFu8NNTuOXIX09HR2795NtWrVzI5SZFWrQt269lHrV60yO42IiHgSFfIiIlJqfWIf446bboLKlc3NIsXz5JNPsmTJEhITE1m+fDm33347VquV/v37mx2tWHLPyqt7vYiIFIcKeRERKZVs2TaqTfoX7VjN/fcZZseRYjpw4AD9+/enYcOG9OvXj0qVKrFy5Uoqe9gRGV0nLyIiV0LXyIuISKm0cfLvjEibwP1Mwu+GZCDQ7EhSDHPnzjU7glPkFvIrV0JODlit5uYRERHPoDPyIiJSKp38j71f/cb6dxFYUUW8mKNZMyhbFtLTYfNms9OIiIinUCEvIiKlzum/T9Ni51cAlHtUc8eLeaxW6NjRfnvFCnOziIiI51AhLyIipc6Gl/+PUNI4YK1JixHXmh1HSjldJy8iIsVleiG/dOlSevfuTUREBBaLhfnz5+d53DAMxo4dS7Vq1QgMDKR79+4kJCRc9nknT55MVFQUAQEBdOjQgdWrV7voHYiIiKexfv4pALvaD8DH1/SmUEo5jVwvIiLFZfq3l4yMDKKjo5k8eXK+j7/xxhv8+9//ZurUqaxatYrg4GB69OjBmTNnCnzOL774gtGjR/Piiy+yfv16oqOj6dGjB0eOHHHV2xAREQ+Rsu0orQ//BEDkc+pWL+bL7Vq/ezccPmxuFhER8QymF/I9e/bklVde4fbbb7/kMcMwmDhxIs8//zy33XYbLVq04OOPP+bgwYOXnLm/0DvvvMOwYcN44IEHaNKkCVOnTiUoKIgZM2a48J2IiIgn+PWDBI5Qha1Bbah7S2Oz44hQvjw0bWq/revkRUSkKEwv5Auzd+9ekpOT6d69u2NdaGgoHTp0YEUBLV1WVhbr1q3Ls4+Pjw/du3cvcJ/MzEzS0tLyLCIi4p3e/COWmiSx4un5ZkcRcdB18iIiUhxuXcgnJycDULVq1Tzrq1at6njsYikpKeTk5BRrn7i4OEJDQx1LZGSkE9KLiIi72bED1qwBi9VK70dqmB1HxCG3kNcZeRERKQq3LuRLypgxY0hNTXUs+/fvNzuSiIi4wA//3o2VbHr0gCpVzE4jcl5uIb9mDWRlmZtFRETcn1sX8uHh4QAcvmjkl8OHDzseu1hYWBhWq7VY+/j7+xMSEpJnERER72LLttFvWjf2E8lj18SbHUckj/r1oVIlyMyEDRvMTiMiIu7OrQv52rVrEx4ezuLFix3r0tLSWLVqFTG5c7VcxM/PjzZt2uTZx2azsXjx4gL3ERER77dpyjJq5OwjmAy6PNTQ7DgieVgsmoZORESKzvRCPj09nfj4eOLj4wH7AHfx8fEkJSVhsVgYNWoUr7zyCt999x2bNm1i4MCBRERE0KdPH8dzdOvWjUmTJjnujx49mg8//JDZs2ezbds2HnnkETIyMnjggQdK+N2JiIi7SJ1snzv+z/p3EVgx0OQ0IpfSgHciIlJUvmYHWLt2Ldddd53j/ujRowEYNGgQs2bN4umnnyYjI4OHHnqIEydO0LlzZxYsWEBAQIBjn927d5OSkuK4f/fdd3P06FHGjh1LcnIyLVu2ZMGCBZcMgCciIqXDmRNniN7xJQBlH9bc8eKeLizkDcN+ll5ERCQ/FsMwDLNDuJu0tDRCQ0NJTU3V9fIiIl5g+T/nEftOX/6yRlLtTCI+vqZ3SCs2tU3O526f6alTEBICOTmQmAi1apmdSERESlJx2iXP+yYjIiJSTNbPPgEgod0AjyzipXQICoJWrey3NQ2diIgURt9mRETEq6Xs/JvWh38EoMYYdasX96br5EVEpChUyIuIiFeb+3MFOrKSf0e8Rr1bm5gdR6RQGrleRESKQoW8iIh4tU8+tbCeNtieesbsKCKXlXtGPj4eMjJMjSIiIm5MhbyIiHitHTtg9WqwWqF/f7PTiFxeZCRUr24f8G7NGrPTiIiIu1IhLyIiXivx0TeYwQM83GEDmoHUu7322mtYLBZGjRpldpSrYrHoOnkREbk8FfIiIuKVbDkGjZdO5QFmcW+rbWbHERdas2YN06ZNo0WLFmZHcYrcQl4j14uISEFUyIuIiFfaNPUPambv5SRlaTmuj9lxxEXS09MZMGAAH374IRUqVDA7jlNceEbeMMzNIiIi7kmFvIiIeKXUSfa54+Pr3UVQWJDJacRVhg8fTq9evejevftlt83MzCQtLS3P4o5atoSAAPj7b9i50+w0IiLijlTIi4iI1zlz4gwttn8JQNmHNXe8t5o7dy7r168nLi6uSNvHxcURGhrqWCIjI12c8Mr4+UHbtvbbuk5eRETyo0JeRES8zoZXfqA8JzhorUH0yK5mxxEX2L9/PyNHjuSzzz4jICCgSPuMGTOG1NRUx7J//34Xp7xyGvBOREQK42t2ABEREWezfGbvVr+z7QAifHXM2hutW7eOI0eO0Lp1a8e6nJwcli5dyqRJk8jMzMRqtebZx9/fH39//5KOekVUyIuISGFUyIuIiFdJSYH/O9KBMDZTY4y61Xurbt26sWnTpjzrHnjgARo1asQzzzxzSRHvaWJi7D+3boUTJ6B8eTPTiIiIu1EhLyIiXuWLL2CCbQw/tXyW9bdZzI4jLlKuXDmaNWuWZ11wcDCVKlW6ZL0nqlIF6tWDXbtg5Uq46SazE4mIiDtRf0MREfEqn9h71XP/QBXx4tlyz8qre72IiFxMZ+RFRMRr7Pl1H9VXrSXAcgv9+3vGtdDiPL/99pvZEZwqNtZ+YEqFvIiIXExn5EVExGvse3EG/+UuFla5l/Bws9OIXJ3cAe9WrYLsbHOziIiIe1EhLyIiXsGwGdRZ8SkAPnfeYXIakavXtCmUKwfp6bB5s9lpRETEnTi1kF+1apUzn05ERKTINk5dTq3sPaQTTMtxfcyOI4U4ffo0f/311yXrt2zZYkIa92W1QseO9tsrVpibRURE3ItTC/m+ffs68+lERESKLHWSfZS7+Lp3ElQ52OQ0UpB58+ZRv359evXqRYsWLfKcBLj/fk0XeDHNJy8iIvkp9mB3/fr1y3e9YRj8/fffVx1IRESkuDLTMmm+/UsAgh5SMejOXnnlFdatW0fVqlVZt24dgwYN4rnnnuPee+/FMAyz47kdjVwvIiL5KXYhv2jRIj755BPKli2bZ71hGCxdutRpwURERIpqwys/0NE4TrJPBNGjrjM7jhTi7NmzVK1aFYA2bdqwdOlSbr/9dnbt2oXFoikDL9ahA1gssGcPJCejQRxFRAS4gkK+a9eulCtXjmuvvfaSx1q0aOGUUCIiIsWR/I39AuLtbQcQ7mc1OY0UpkqVKmzcuNHxnaFixYosXLiQQYMGsXHjRpPTuZ/y5e2D3m3ebL9O/vbbzU4kIiLuoMjXyJ88eRKAr7/+Ot8iHmDhwoXOSSUiIlJEx45Bv31v0ohthL/6uNlxpAC53yM++eQTqlSpkucxPz8/Pv/8c5YsWWJGNLen6+RFRORiRS7kr7nmGpKTk12ZJV9RUVFYLJZLluHDh+e7/axZsy7ZNiAgoIRTi4hISfnySzh7FgJbNqJR9xpmx5EC5H6PqFGjBuEF9A/v1KlTCafyDLmFvEauFxGRXEUu5Fu1akWHDh3Yvn17nvXx8fHcfPPNTg+Wa82aNRw6dMix5J71L2yE/JCQkDz77Nu3z2X5RETEXPNmpQOgAc/dm1nfI7xBbiG/di1kZpqbRURE3EORC/mZM2cyePBgOnfuzLJly9i5cyf9+vWjTZs2WK2uux6xcuXKhIeHO5bvv/+eunXr0qVLlwL3sVgsefbJHVRHRES8y77Fu/h+dWXmcg/977aZHUcKYdb3CG9Qrx6EhdmL+A0bzE4jIiLuoFiD3Y0fPx5/f39uuOEGcnJy6NatGytWrKB9+/auypdHVlYWn376KaNHjy50ZNv09HRq1aqFzWajdevWTJgwgaZNmxa4fWZmJpkXHOJOS0tzam4REXGNvS9/Si3OUK/ScapVL/KxaTGJ2d8jPJXFYp+G7v/+z36dfMeOZicSERGzFflbz+HDhxk5ciSvvPIKTZo0oUyZMgwePLhEG9/58+dz4sQJBg8eXOA2DRs2ZMaMGXz77bd8+umn2Gw2YmNjOXDgQIH7xMXFERoa6lgiIyNdkF5ERJzJsBnU/uNTAM70Vb96d+cO3yM8mQa8ExGRCxX5jHzt2rVp2LAhX331Fb169WLBggXcfffdJCUl8dRTT7kyo8P06dPp2bMnERERBW4TExNDTEyM435sbCyNGzdm2rRpvPzyy/nuM2bMGEaPHu24n5aWpmJeRMRkSUmQklLw44f+u4Je2btJJxifO29n/Xr7+rAwqFmzZDJK0bnD9whPdmEhbxj2s/QiIlJ6FbmQnzFjBvfcc4/j/k033cSvv/7KLbfcQmJiIpMnT3ZJwFz79u1j0aJFfP3118Xar0yZMrRq1Ypdu3YVuI2/vz/+/v5XG1FERJwkKQkaNoQzZwre5j98AsDX3MGgG4Id6wMCYMcOFfPuxuzvEZ6ubVvw9YVDh+z/P2rVMjuRiIiYqchd6y9sfHO1bt2a5cuX88svvzg1VH5mzpxJlSpV6NWrV7H2y8nJYdOmTVSrVs1FyURExNlSUgov4v3I5G6+AOAT8narP3Om8DP5Yg6zv0d4uqAgaNXKflvd60VE5KpHBoqKimK5i1sUm83GzJkzGTRoEL6+eTsRDBw4kDFjxjjuv/TSS/zvf/9jz549rF+/nvvuu499+/bx4IMPujSjiIiUnJv5kYoc5y8i+IXrzY4jV6Ekvkd4i9wrB/VxiYhIsUatL0iFChWc8TQFWrRoEUlJSQwZMuSSx5KSkvDxOX884vjx4wwbNozk5GQqVKhAmzZtWL58OU2aNHFpRhERcb16JDCEGTRiO4u5nmV0woamLvN0rv4e4S1iY+Hf/1YhLyIiYDEMwzA7hLtJS0sjNDSU1NRUQkJCzI4jIlLqrF8PbdrkXTeYmXzEgxhYsGA4fg5lOrMZnGfbdeugdeuSy1sS1DY5n6d9pvv328d+sFrhxAkoW9bsRCIi4kzFaZc06a6IiLi9eiTwEQ9ixYYvOY6fPtiYzlDqUvCApuKdpkyZQosWLQgJCSEkJISYmBh++ukns2O5VGQk1KgBOTmwZo3ZaURExEwq5EVExO0NYQYGl863ZQEMLAxlesmHElPVqFGD1157jXXr1rF27Vquv/56brvtNrZs2WJ2NJfKnYZuxQpzc4iIiLlUyIuIiNuLIhELBV0JZhBFYknGETfQu3dvbr75ZurXr0+DBg149dVXKVu2LCtXrjQ7mktdOJ+8iIiUXk4Z7E5ERMSVEokq5FHLZR4Xb5eTk8NXX31FRkYGMblDu+cjMzOTzMxMx/20tLSSiOdUF56Rt9nAR6dkRERKJf35FxERtzefPvhgu2S9AVgwmM7Qkg8lptu0aRNly5bF39+fhx9+mG+++abQWWri4uIIDQ11LJGRkSWY1jmioyEgAP7+G3buNDuNiIiYRYW8iIi4OYN/8eq56+EhGyvZ+JCNFRs+DGU6u6lndkgxQcOGDYmPj2fVqlU88sgjDBo0iK1btxa4/ZgxY0hNTXUs+/fvL8G0zuHnB+3a2W+re72ISOmlrvUiIuLWHmYqt/J/ZOLHncyjE8uJIpFEopjOUBXxpZifnx/16tn//du0acOaNWt47733mDZtWr7b+/v74+/vX5IRXSI2Fn7/3V7IDxlidhoRETGDCnkREXFrddkNwLO8xg/05gd6m5xI3JXNZstzDby30sj1IiKiQl5ERNzaU7zFfPqwnFizo4gbGTNmDD179qRmzZqcPHmSOXPm8Ntvv/Hzzz+bHc3lcsfz27oVjh+HChXMzSMiIiVP18iLiIjbCQsDq8/56eb+oDNGEZusgAD7/uLdjhw5wsCBA2nYsCHdunVjzZo1/Pzzz9xwww1mR3O5ypWhfn37bS+fbU9ERAqgM/IiIuJ2jn72PxbY3mAws3jy3Rpce23R9w0Lg5o1XZdN3MP06dPNjmCqmBhISLBfJ9+zp9lpRESkpKmQFxERt5KyPYUazw+iDcnMjn6XbqPeNjuSiNuJjYWPP9bI9SIipZW61ouIiNswbAa7uw6lqi2Z3X6NiVn0stmRRNxS7oB3q1ZBdra5WUREpOSpkBcREbexdMA0Ohz+jkz8yP7kc4LCgsyOJOKWmjSBkBDIyIDNm81OIyIiJU2FvIiIuIXd32+j3dzRAKy47TUa9os2OZGI+7JaoWNH+211rxcRKX1UyIuIiOky0zLJ7tufIE6zttKNXDtvpNmRRNxebvd6FfIiIqWPCnkRETHdm/9MJudMFimWMCIXzcLHV82TyOXkzievQl5EpPTRNyURETHV//4HL3xUi7asZfM7/6Nqy2pmRxLxCB06gMUCe/fCoUNmpxERkZKkQl5EREyTctRg0CD77cGPBNF1VCtzA4l4kNBQaNbMfnvFCnOziIhIyVIhLyIipjBsBjuj+zIo+TWaNsrhrbfMTiTieXKvk1chLyJSuqiQFxERU/x+3zRiD/2X8bzIVxMSCNJMcyLFpgHvRERKJxXyIiJS4nZ/v422n5+bau7WOBrf3sjkRCKeKbeQX7sWMjPNzSIiIiVHhbyIiJSozLRMzva79/xUc/8dZXYkEY9Vty6EhUFWFqxfb3YaEREpKSrkRUSkRK24/l80Oh3PMUslTTUncpUsFnWvFxEpjdz+29O4ceOwWCx5lkaNCu+C+dVXX9GoUSMCAgJo3rw5P/74YwmlFRGRwqx7bSFd170NwK4xMzTVnIgTqJAXESl93L6QB2jatCmHDh1yLMuWLStw2+XLl9O/f3+GDh3Khg0b6NOnD3369GHz5s0lmFhERC6WkgKfvbafTPxY2uRhOrx6q9mRRLzChYW8YZibRURESoav2QGKwtfXl/Dw8CJt+95773HTTTfx1FNPAfDyyy+zcOFCJk2axNSpU10ZU0RECmAY8OCD8G3qEHbXbsfnS+qaHUnEa7RtC76+kJwM+/ZBVJTZiURExNU84ox8QkICERER1KlThwEDBpCUlFTgtitWrKB79+551vXo0YMVhUywmpmZSVpaWp5FRESc54NpBt9+C35+MP7r5gSFaa45EWcJDITWre231b1eRKR0cPtCvkOHDsyaNYsFCxYwZcoU9u7dyzXXXMPJkyfz3T45OZmqVavmWVe1alWSk5MLfI24uDhCQ0MdS2RkpFPfg4hIabb7h+3EPNqKNqwlLg5atjQ7kYj3iYmx/1QhLyJSOrh9Id+zZ0/69u1LixYt6NGjBz/++CMnTpzgyy+/dNprjBkzhtTUVMeyf/9+pz23iEhplpmWSVbfe2lh/MnksBcZNcrsRCLeSQPeiYiULm5fyF+sfPnyNGjQgF27duX7eHh4OIcPH86z7vDhw4VeY+/v709ISEieRURErt6Kbs/T+PQGjlkqUevnD/HxuFZH3FVcXBzt2rWjXLlyVKlShT59+rBjxw6zY5kmt5D/809ITzc3i4iIuJ7HfaVKT09n9+7dVKuW/5RFMTExLF68OM+6hQsXEpPb50xERErE+jcW0XXtWwDsHjOd8NYRJicSb7JkyRKGDx/OypUrWbhwIWfPnuXGG28kIyPD7GimqFEDIiPBZoM1a8xOIyIirub2hfyTTz7JkiVLSExMZPny5dx+++1YrVb69+8PwMCBAxkzZoxj+5EjR7JgwQLefvtttm/fzrhx41i7di0jRoww6y2IiJQ6x3Yeo9qYQQAsbfIP2r96m8mJxNssWLCAwYMH07RpU6Kjo5k1axZJSUmsW7fO7GimUfd6EZHSw+0L+QMHDtC/f38aNmxIv379qFSpEitXrqRy5coAJCUlcejQIcf2sbGxzJkzhw8++IDo6GjmzZvH/PnzadasmVlvQUSkVDFsBgldHqSa7SC7/RrRdsk7ZkeSUiA1NRWAihUrFriNt89So0JeRKT0sBiGYZgdwt2kpaURGhpKampqsa6XT0qClJTiv15YGNSsWfz9RETc0Yz3M6j2+F10YzF75qyiUf9WZkfyClfaNpUGNpuNW2+9lRMnTrBs2bICtxs3bhzjx4+/ZL23fKZr1kD79lChgv37iMakEBHxLMVp61XI5+NKviwlJUHDhnDmTPFfLyAAduxQMS8inm/7dvt81mdO25g9Kp77321tdiSvoUK+YI888gg//fQTy5Yto0aNGgVul5mZSWZmpuN+WloakZGRXvOZnj0LoaFw+jRs3QqNG5udSEREiqM4bb2O1TpJSsqVFfFg3+9KzuSLiLiTzNM27r3XXkR06+7DgLdVxIvrjRgxgu+//55ff/210CIevH+WmjJloF07+211rxcR8W4q5EVExClWdHmWERuGULNiOrNnq1uvuJZhGIwYMYJvvvmGX375hdq1a5sdyS3kXie/YoW5OURExLV8zQ4gIiKeb/2bi+m65k26Ao0fvYuIiJvNjiRebvjw4cyZM4dvv/2WcuXKkZycDEBoaCiBgYEmpzOPBrwTESkddL5ERESuyt8Jx6j27EAAljb+BzEvq4gX15syZQqpqal07dqVatWqOZYvvvjC7Gimiomx/9y2Df7+29wsIiLiOjojLyIiV8ywGezs8iAdc6eaW6qp5qRkaKze/IWFQf36kJAAK1fCzTquJiLilXRGXkRErtiyQR/S8dB8sijD2VlzCAoLMjuSSKmn7vUiIt5PhbyIiFyRvT9tp82nowBY3muC5osXcRMq5EVEvJ8KeRERKbasLIgbmUw6ZVlXsTvXzh9tdiQROSe3kF+9GrKzzc0iIiKuoUJeRESK7fnn4cOErnQpv5Hqiz7Gx1fNiYi7aNIEQkIgIwM2bTI7jYiIuIK+eYmISLEsXmjjzTftt+NmhhPeqpq5gUQkDx+f86PXq3u9iIh3UiHvJj7+GHJyzE4hIlK4vxOOEXFzNHfxFQ89BH36mJ1IRPKj6+RFRLybCnk38d579qPnmzebnUREJH/2qeaG0Th7M6/7vcA7r2WZHUlECqAz8iIi3k2FvJsIDoY1a6B1axg/3j6QlIiIO/l98Ed0PPQNWZQha9bnBFfwMzuSiBSgQwewWCAxEQ4eNDuNiIg4mwp5JwkLg4CAK9s3IAB++QV694azZ2HcOGjb1l7Yi4i4gz0/7aDNJ6MATTUn4glCQqB5c/vtFSvMzSIiIs7na3YAb1GzJuzYASkpxd83LMy+/7ffwhdfwGOP2UeZ7dgR/vlP+xn6wEDnZxYRKYqs9CzO3DWAYE6xvkI3TTUn4iFiY2HjRnshf+edZqcRERFnUiHvRDVr2pcrZbHAPfdAt24wciR8/jm8+SbMnw8ffQTXXuu0qCIiRba82wt0PbWOvy0Vqfa/2ZpqTsRDxMbC1Km6Tl5ExBvp25gbqlwZ5syB776DiAhISIAuXWD4cDh50ux0IlKa/LLYYOPqMwAkPP0R1dpWNzmRiBRV7sj169bBmTPmZhEREedSIe/GeveGLVvgwQft9//zH2jWDH7+2dxcIlI6HDsGAwdZGMl7jL9rIx1eu93sSCJSDHXq2E8OZGXB+vVmpxEREWdS13o3V748fPihvcv9sGGwdy/cdBMMGgTvvAMVK5qdUEQ8TVLS5cfzMGwGzzxl46+/rNSqBd1HNScp6eouHxKRkmWx2M/Kf/utvXt97hl6ERHxfCrkPUS3bvYB8P71L/j3v2H2bFiwACZP1gA2IlJ0SUnQsOHlu9kOZTrP8yk7+IR9+yLp3Nk+w8aOHSrmRTzJhYW8iIh4D3Wt9yDBwTBxIixbBo0aweHDcNdd9iU52ex0IuIJUlIuX8Q3YAfvMZKuLOEu5jnWnzlzZTNziIh5cs/Cr1gBhmFuFhERcR4V8h4oNhY2bLCfnbda4b//hSZN4OOP1UiLyNUpQxafYZ9qbhHdeI+RZkcSkavQpg2UKWM/4J+YaHYaERFxFhXyHiogAF55BdasgZYt4fhx+3XzN99s7zorIlIc9UhgAmOIpyVtWcdxQhnEbAw1EyIeLTAQWre231b3ehER7+H239Di4uJo164d5cqVo0qVKvTp04cdO3YUus+sWbOwWCx5loCAgBJKXLJatYLVq2HCBPD3t18337QpTJkCNpvZ6UTEEwxmJttpxFO8QWO2ARBKGjew0ORkIuIMMTH2nyrkRUS8h9sX8kuWLGH48OGsXLmShQsXcvbsWW688UYyMjIK3S8kJIRDhw45ln379pVQ4pJXpgyMGQPx8fZu9+np8OijcN119jnoRUQKUo8EPuJBrNjwxYbl3HoLBtMZSl12mZpPpDBLly6ld+/eREREYLFYmD9/vtmR3FLudfIq5EVEvIfbF/ILFixg8ODBNG3alOjoaGbNmkVSUhLr1q0rdD+LxUJ4eLhjqVq1agklNk+jRrB0Kbz3HgQF2W+3aAFvvQXZ2WanExF3NIQZGI7y/TwLYGBhKNNLPpRIEWVkZBAdHc3kyZPNjuLWcs/Ib9wIJ0+am0VERJzD46afS01NBaDiZSZQT09Pp1atWthsNlq3bs2ECRNo2rRpvttmZmaSmZnpuJ+Wlua8wCXMaoXHH4feveGhh2DRInjqKfjyS5g+HZo3L9oc0vkJC9O0UyLeoirJjOYd6rIbCwWNkmkQRWJJxhIplp49e9KzZ0+zY7i9GjXs7XdSkn1sneuvNzuRiIhcLY8q5G02G6NGjaJTp040a9aswO0aNmzIjBkzaNGiBampqbz11lvExsayZcsWatSoccn2cXFxjB8/3pXRS1zt2vC//8GMGfDPf9ob7jZtYPhw+/XzFxy3KDLNIS3i+dJ2JvM2b/AIUwjkDL/TKd8z8nYWEokqyXgiLuVNB+6LKzbWXsgvX65CXkTEG7h91/oLDR8+nM2bNzN37txCt4uJiWHgwIG0bNmSLl268PXXX1O5cmWmTZuW7/ZjxowhNTXVsezfv98V8UucxQJDh8LWrXDrrXD2rH0e+isp4kFzSIt4siN/HmJJ6yfo0L82o3mXQM6wgo58yDAsGJeckzc4f528iLeIi4sjNDTUsURGRpodqcToOnkREe/iMWfkR4wYwffff8/SpUvzPatemDJlytCqVSt27cp/0CZ/f3/8/f2dEdMtRUTA/Pn27vUPPwwnTpidSERKyqGDBttufZqYdZPowhkAVtCRFxnPQm4ALPicK9htWM51s7f/HMp0dlPP1PwizjRmzBhGjx7tuJ+WlubVxfyFl9JVqGD/+fvvsHYt+BRyKkeX0omIuD+3L+QNw+Cxxx7jm2++4bfffqN27drFfo6cnBw2bdrEzTff7IKEnsFigbvvtjfO3btf+ng9EhjCDKJIJJEoZjCEXdQv+aAi4hSHDsEbb8DUqRYmnznG9Zxhc9mO7Bk0ntsm2wv4XLMZzDI6M5Tpjr8B0xmqIl68jrcfuL9QUhI0bGjvTXeh9HRo167wfXUpnYiI+3P7Qn748OHMmTOHb7/9lnLlypGcnAxAaGgogYGBAAwcOJDq1asTFxcHwEsvvUTHjh2pV68eJ06c4M0332Tfvn08+OCDpr0Pd5F7RP5Cg5nJRzyIce4snIGFp3mDoUxnNoNLPKOIXLnD8YfY9sAbjN76IBuy7AN8/tDqBaL73kPrZ24gK94C+QzwvZt6PEdcCacVEVdJSbm0iC+q3EvpVMiLiLgvty/kp0yZAkDXrl3zrJ85cyaDBw8GICkpCZ8L+ogdP36cYcOGkZycTIUKFWjTpg3Lly+nSZMmJRXbY1w4h/SFDGA6Q1lGZ52VE/EAh+MPsf2B12kfP42unOFpDvF+7FzGjYPu3WtjsRS/N5OIu0tPT89z2dzevXuJj4+nYsWK1FQVKiIiXsztC3nDKGhapPN+++23PPffffdd3n33XRcl8i6Xn0P6I57jtZIPJiJFkrz+IDuGvE6HP6fRBftIlhvLxdLgX0NZ9rT9spoLhYXZu81eyZm6gAD7/iLuYu3atVx33XWO+7nXvw8aNIhZs2aZlMq96VI6ERHv4PaFvLhWFIkFziHtSw6P8T7+ZDGPu1hJRzxsogMRr3XwIPx5x3iuWxWXp4DPfn48rZ7shsUn/ynlata0X/t6JTNQaAAscTddu3Yt0gF/sdOldCIi3kOFfCmXSFSBc0gbQFlOMZp3Gc27rKM1r7++jkcegWuuAau1ZLOKCPz1F7z+OnzwATyeGUhPMtlYrhPZz48rtIC/UM2aKshFShtdSici4l10erWUm8GQAueQtuHDw/yHTxlAKiGspj1ffgnXXQcR4TYWtRjNurj/cfbUWTOii5Qqh9b+xZLox3ks6v94/33IzIQNHR9l/esLaX7id1o/3b1IRbyIlE6Xv5RuesmHEhGRK6Yz8qXcLuozlOkFziE9m8FM4xH8yKQs6fTuDX/8AfVTltM95V3Y9C7H/1WBzXVvw//eu4j+Z3f8Q0rH1D4iJeHQ2r/YOfQ1Omz8kC5kEsLvHO10C+PGW7j++rJYLPnMJykicpHCLqWzYBBFYskGEhGRq6Iz8sJsBtOQHbzFU3xFP97kKRqyI8/1cln48zeVGDcOkpPhjY8qsbTxPzhqqUIF4zjX7JpF+5du4UxoFf6IGsDidzdy6pRpb0nE4x1a+xdLWjxGxXZ16LJxEgFk8mdIZ3jrbZYuhW7dLh3ITkSkIIVdSmfFRieWcRvz8SGnhJOJiMiVsBgaJeYSaWlphIaGkpqaSkhIiNlxnGr9emjT5sr3X7cOWrc+fz8nK4dNU5aROn0eDbZ8TTXbQQBi+YM/g2Lp1QsGdEum221lKRte9irTi7ivpCTnDCD311+w6p536bXsWfzJAuDPkM7Yxo6n5RPXqft8KebNbZNZvPkzvbi9r0cC22mED7Y85Xzul8DcdZtoRkvi+W2plWuuKaGwIiICFK9dUtd6uSpWPystR3aBkV2wZb/HphmrSJ7xI4cOdeRUEnz1FXT76kWsD3/MqvAenL3tLpqN6U35WqFFen5nFUcirpSUBA0bXvmUbjt2gI8PvPYafPgh9Miqyx1kqYAXEacp7FK6J3mLyhzlH0zjDzphw8rNN8PDD8PIu5Op0Tbc7PgiInIRnZHPhzcfoXdGwVGUAtkw7Gfv//tfuHNiZ9qe+cPxWBZliK98A5m97qTpc7dRsX4lU7OKXK2r6elSnQN8VO81Fu+ty1s5TwBw7TUG7/ZdTqvhsSrgxcGb2yazePNnWtDfpbrsYijTHfPIT2eoY7T6QE4RxCmOEQZAW9awghhWR95JyIujaTa0Q0m+BRGRUqc47ZIK+Xx4c8MOJX+W27AZJHy9ib/+/V8iV8+jXuZWx2Pbachj3bdz551w++1Qter5/Zx9GYCIq1zJ72p1DvAsrzGMD/EnixQqcW+nJMa8HETXrrr+XS7l7W2TGbz5M73aNvTdd6Hc+68ydM/zjnWbysaQPuwJ2k24Hd8AdeoUEXE2da2XQpX0HNIWHwsN7mpBg7taAOPZ/f029k/8L+F//Jdvz9zIokWwaBGMeiST3yrcTuZ1N9Hg2TvAWqPkQoq4QD0SGMIMx5mvGQzhDAE8y2s8yEeOa+DXlb0W60vj+N+oQAoYi0pEpERdey20HvUvdnx5C0efe5f2u+fQPH0FvLuCA/+uxa6bH6fVh48SWjXA7KgiIqWSzsjnw5uP0LubXduz+fo7X+bNg7A1P/IjvRyP/RnYkY9P38V/uZN9RDnW51cc7aL+Jc+tM/JSUvI78zWYmXzEgxjnrkE1sOCDjRx8KHNuVOglXMuLjOeddV31uyqXpbbJ+bz5M3X25WlHNiazdcR/aLZsCmFGCklEEh28m0EPluHxx6FOHedlFxEprdS1/ip5c8Puzg6sOcSuV+dS8Zd5tDi5PM9ja2nDKCZSn4RLiqML57y/kAp5KSkFjQ5txZZnu9w/tqtpzzO8zhK6AvpdlaJR2+R83v6ZuuJSutN/n2btqE/5YXEArx+8HwA/y1m+r/UoVZ8aRPOHO2lsDxGRK6Su9eKRarSrRo35TwBPkLz+ICuf/YbQhfO4lqW0ZR0hpPIRD+ZbHE1nKMvo7BiwR8QMVrJpxmZe45lzI0LnZQGysfIL1zuKeBERV3HFpXSBFQO55uNhdDbguv/BO+9Ahf/9lxsSP4LhH7HlqXakDnmCdq/fRZmgMs59cRERcVAhL24pvHUENV8bTpuFw6nMEW5gIdfwO0Y+FxBbAAsGn3EvH/APttKE1bQHrCWaWVPllT7HjsHKlXBq3Bss4mc6sIqyZFxmL4MoEksinoiIy1gs0KOHfUn4vhVLn3yQ9js+oempNTDpXg5OeZodNz5Gq8nDKF+7gtlxRUS8jgp5cXtHqcIcBjCH/vme5QTwwaADa+jAGs7gTzAZ9OgB0dEw2DKbyLDThHZsTI0bmxDWuLLTM2qqPNdwp4Mjtmwbu7/fxqGvl3Nq/Q5GZb/Fjh32x37iF7rxCwCphHCMikSxD598f18tJF4w5oOIiKerf0tD6t/yIUe3vMrKx6bS9LfJROQcIOKnZ0iv8xL/GrSDwf+qTv1Lh7MB3OtvvYiIp1AhLx4jkah8z8gD5ODDelpznArY8MGGlZQUWLwY3mQirYiHufZtUyxhHCjXhNQaTbA1bUHOQ4/QuDFERFz5lF8pKVdWxIN9v5QUfRm5mNkHR1L3p5Hw6SrS/7ecsptWUP/YSuqT6hhWMZV/AtVo0ADWhjzC12vvYAUxbKUJddjDdhqdG8XhPAN775HpDL3yYCIibqpy0yp0/WUsmWnPsGz051T57B0OnwllwuzqxH0MvXvDc30TaH9vPcd19Gb/rRcR8VQq5MVjzGAIT/PGuWLovNxznv35PM818rNmgWHAsQ96s3p3dcL/3kqN7ETCjBTC0pbC1qVs29qIJl89AkBICHzi/yAVyhvkNGhMcLsmVOvWhIiONfHx9Smpt+lynnLmoyQPjhg2g70/7+T3A7X5Y40fK1bAI5vH8Cj/ybNdBkHsKN+B1CYxfPywhVY97Z/L+vW35Rnsbhf1Gcp0pjMU27kBGblgYEaN5SAi3sw/xJ/OHw3G+GAQh+b/Ta8Z8MMP8Md3KbT4rgU7HmpEysDRtH/rblJS/DzqQLintKHgOVk9JScoq6t4SlZ3y6lCXjxGcYuj5s3PjQQ++CXHulMpp0hauINjv2/l7J9b2Z8aQsNs2LUL0tIMuvIlIUdPQgLwAzDOXrwlBTZiT9T1bLz/TRo3hiZNoE6UDV+/wgv8ok6VV1J05sMuPTmdnZ+tIW3BcoI2rqDe0RXUMf7mPv5gBbEALCeGXr4LOFAjhux2MVTpE0v9O5rTOqBofzZnM5hldGYo0x3//tMZqiJeREoNi4+FLndUossdsH07LPrnGowfLTQ6HQ/TBpL84TMcaDucijzM31QyO+5leVIb6ilZPSUnKKureEpWd8ypQl48ytUWR0FhQTTq3wr6t3Ksux/IzISE7Ta2zJ5O5oat+O3eSuUjW6mZuZNgTtH49Hp2b6vGc8/l7mVwgEgyAiqRUrkJx6o05i6asJUmJFCfs/jlO4/407yR71R5JcUbLgEo7sERw4A9e2DFCvj7q8XcsPApGpz+k9YXzX5wBn9uarCXzrfFEhMDMR0HEF7tPmpdRdbd1OM54q7iGUREvEOjRtDoh578nbCf1SOm0WjRJMJth7h19fPs51U+ZiAvMZZDRADudyAcPKsN9ZSsnpITlNVVPCWrO+ZUIS8exxXFkb8/NIu2wjt986zPPpPN3iV7OPLbVtKOhnBfJmzdCilbj1L9zEE4c5AG+zfBfuiduw9W5nMbtzNfU+U5WVEOjgRwmras5cBjyzm7ZwVTTg9idurtAMQQyONsAOAvayT7ImI42zaWsN4x1O/bkrFl/S54Nc2DLCLibBXrV6Lrz8+Rlf4kfzz1JRVmvkOTzA08wExeZDxQtL/1IiKlnQp5kUL4BvhSu0cDavdoQAfg3nPrbdlhHFi5l0OLt5KxZivZG7cSvH8bTdhKKGnUIqnAqfJ8sLGRFiQSxVEqM58+vPrqEzRpAlUq5dBm95cERFYmuFYYofUqU7FhZfzyFJjO545nPi5WjwQ+4sF8D47MYAhd+ZUmbKMVGyhDNiy3P76W6swpczutW0Pn9q1ZXuYrou6JoXq76lR3UrawMHu3qSvtbhUW5qQgIqXQ5MmTefPNN0lOTiY6Opr333+f9u3bmx1LLsOvrB+dptzHuqED6NJuKc3ZxBGqFvq33l0PhHtCG5rLU7J6Sk5QVlfxlKxm5lQhL27LnYsjH18fanSOokbnKOBm1q/n3GBnBhEcZDLDac36fPe1AEGcpgnbgG1spAUTv4avv4bKHOOI43DBeamEcMK3Mour3sv81i9RuTJUrZTNjdvewzc8DP8alQmOqkxI3cpUaFCZ4CrBRX4vZp/5yMnK4fSxU5zKMMiwhpCRAacyDA5/uYSbySCIUwSTQX8+z3d/C/aR4AfzsWPdQaqxq3IMtvYxdLqrO2n32H8nIAC4y+nvoWZN+7VP7jQAikhp8MUXXzB69GimTp1Khw4dmDhxIj169GDHjh1UqVLF7HhSBBYfC0vpwlK6ADCEGYUeCF9Ed8byEp8wEIBprxyl/+GJULYslrLBWELK4hsSjG/5spQpH4xP3dqUqVeL4GAoG2QjuEwW/iH+jlHzr5bZbWhxeEpWT8kJyuoqnpLV7Jwq5MVteWZxZOEg1dlGY27h+3y3yMbKLAbzGQOozFH2UpuHHgJfX8jZl8WGP7pS9vRRymcdpaKRghUboaQRmp1G+l8n+L+/7M9ThWNM4Ml8X+MUgXxVdgj/bjDJUfQ/sP1p0sqE8SCVOXpuCSLjMmc+OrGfmgSfK6hTlp9i89FgTgRX59QpOJOaSaVf55GTfgrjZAZGximMjFP4nM7AcvoUCSFt+KH6Q2RkQHb6GSZtiMU/5xQBORkE2E4RZGQQQCZlgQXcSV/mOTLkcD29852H/VI2IIH6jGM8y4kliZqsW2CxD3ZYQmrWVEEuUtLeeecdhg0bxgMPPADA1KlT+eGHH5gxYwbPPvusyenkSkSReG4w20tZgCj2Ecl+x7qV3xxkGhMKfL7XeZpneR2A2iSyh7pkYyWDYE77BHPaWpZM32Ayy5RlWcTd/NJkBMHBULHMSfpsfAmCgyE4GEu5shzPCuYOgskgmL3UxoZPoW3ocmLYRf1zBybMvVzLU3o6eEpOUFZX8ZSs7pDTYwr54nad++qrr3jhhRdITEykfv36vP7669x8880lmFicwVOLo8KmyrNg8BrP5vnPPe0f50bYpwbwq2O9LdvG8X0nOL7zKOl7j9LYVplpfnD0KJxJNFi2aAABJ48SfOoooVlHqZhzlAAyCeI0J9KtrD/XKaAKx5jNuwDcWoT8uWc+dtKAPOPyPwYfMZRhfARAOTJJ474Cn+cQ/fiGh849px+Nz12fnp9gMggKgqAgCA62sOVQKzKzLGQQzCmCqM1eGrATn3y+4Nmw8jV3Mpf+RXh3IuINsrKyWLduHWPGjHGs8/HxoXv37qxYsSLffTIzM8nMzHTcT0tLc3lOKZ5EovI9Iw+Qgw/fcwvfcptjXfe7KrBk52P4nM7AeiYd38wMymRm4H82Hf/sDE4FVaeiD2RkQHBmBgC+5NgPkNvS7EeCzwKn4X9pHfl6u/15a3GMibx1SYY+535O5R8cp0KBvQes2NhJIwBmMpghzKRNGyhLOscpj4El32W+9S4eCvgYiwV8LAZJJytgWPJuY8MHsLDE/wYeK/8pFgtYLLDscH38yMTAQmXDwp4L9jlNQKE9HVbSgaPYe7FspAV33vkFgYH2bWbs60bl7EP5/pvs86vPYzW/ddyfnNSbGmf35LvtYd8a/CPqZ8f9t/bfTe0zW9l4wTZVScbnosLowpzLiSWZcMd6vxjYZA1kQL1VjnXP//UILU/9kW8GAwt96//puP/kwdF0zFiU77YA99ZdTZZPAAAPH/gXf/J/jsfCC8kKRr5Zd577UvVw1M+klKkGwMCjb3PridkFZhhZaz5/+dUB4O5jk+n797QCt30m8nN2BzTl9Gl4gBmMYuJlsxpYeJt/Upu9+WYFeDXiP2wI7gzA9anf8OiRFwvM8Fb426wsdwMAnU4u4InkpwvcdlLVV/ipjP2b6bUs4X0eu2zWoUxnPn34kGH5ZgX4OGw031UYDED90xuJO1Dwd9UvKj7KV5UeBqBmZgLvJN2Z73ZBNhjFA0zkCQBG8t5lc7p6wGOPKOSL23Vu+fLl9O/fn7i4OG655RbmzJlDnz59WL9+Pc2aNTPhHUhp46x5xH18fahQtyIV6lYEGl70aDjwaZ41hs3gZHI6x3ce5Zr0AP4Pe4+Gk4m+/PbjU2QnH+XM/qOEkUJljlKTJPv15Pm4+PxBFmU4RRBlgvyoX91ecIcEBbF+RzfO+gWT7R+MzT8IW2AwRmAQRlAwATWbMaWD/YRGUJAPa7cvoExIIGVCg/CvGIx/hSACw4IJrBTETRUDybjgD/H69etod8Hc7PVIYDuNznVcuuA9Yz84Mp2hRfpMRcQ7pKSkkJOTQ9WqVfOsr1q1Ktu3b893n7i4OMaPH18S8eQKFXYgHOCfvJ2nDR0wpiatW/+7wOcbf24ByM5qSmryCU6nZHA6JYPMY+lkHc8g63gGZ0+k0yi0HpMr2ot+40gwv/32JD6n0vE5k4H1TAaWjHSyUzMIJoNEoojmzwJ7D1zo4lbLl5wCt7XmZJJhP96ADzZCSaWglwg4c4Lk5PP3q7GfADLz3fYwVQrt6RDG34TxNwBphJCYeP7xCBKoeUEviAudzvJh69bz96uyi3rk///PJ/M0W7acvx/GbhqyOd9tC8pZhaNU4ej5lVmQTjCbNp1fVY5EGrDpkv0BbFjybBtEUoHbAmzZYpB7hWdZ9tOikG2LkjXXru1nL/hEDxWaIXFHpuMTvYvDhW67P+G049E7OFrEvAbVOZB326y8WxzZc9LxaFuOF5rh78RUx6ONSS102xP7TrDr3O1ynCxCXoMoEgkm4/y2WZdudeZACpsO2G8HcrrQDGcPHmHTQfttG5mFbhvO+f9wtdhXSD8be05XsxiGUbS+qybq0KED7dq1Y9KkSQDYbDYiIyN57LHH8u06d/fdd5ORkcH335/v2tyxY0datmzJ1KlTL/t6aWlphIaGkpqaSkhIiPPeiHit89fI51WXXUWaKm/dOkqkG/jFOScwhqd4M98vFdlYmcIjPM8rnCKIbMqUWM78sgIMYlaBB0cuvhapJLOKlAS1TXkdPHiQ6tWrs3z5cmJiYhzrn376aZYsWcKqVasu2Se/M/KRkZH6TE3kSX/ri9uGTmIEr/A8mfiTTjl+/hmaN7VhOZyMYTPAMDBs5xdsNrL9g8mpVAXDsB+c992z85LtcpecoHJkRUTZtzUgYNsGjOwcMAz27jF46aXz5/GHMp3BzMo3aw4+fE5/PjjXgy6dsjzyQWsaNLA/Xm7LSnzO5j1AkFs92AKCONmonWN9ua2rsGaeyrNNLptfAGlNz/9fLbttDYd2nmTie+e3eZCPuJsvLumunJvzC/oxnQcd60aNhIhIKydadj3/vLviKZN27JL9cx1v3e38tns2UubE0YK3je4KVisAKUu28MHL5wu5oXzE3XyZb9ZsrHxF3zwnGR5/HGqcG2X3RLPO2PzsZ/qDDuwk4EhSgRlONI7BFmgf+yjw4G4Ck/cWuG1qw/bkBIeQkACvP5pIvXNl8lCm06+QrFP5B/O5Pd+sACfrtuRsqH3wKf+UvwhO2lZghvTazcmqYD/A6vd3MmUTCz5Yk16rKVv+rsajj0IYR4nmz8tmfZOneJOnHGNRXZwV4FT1+pypap9A2JqRSuiONQVmOF2tDqer2Xs8WE+nE7ptZb7bHfgL/vnvWo6B7N7gSUbzbqE5LzwjX9S/VcVq6w03l5mZaVitVuObb77Js37gwIHGrbfemu8+kZGRxrvvvptn3dixY40WLVoU6TVTU1MNwEhNTb2SyFIKrVuX24xe2bJunTk567HTyMbHsF0UyAZGNj5GXRJMyZlf1tylLgnGBJ415nCPMYFnL8loRlaRkqC2Ka8r+X5wMX2m5vOkv/We3Ia6a1ZPyamsylpSOYvTLvkUXuabr7Cuc8kX9ie6QHJycrG2z8zMJC0tLc8iUhrkXgJgw4ezWMnGh2ys2PAp1iUAJWk39XiOOO7lc54jzi0ziojr+fn50aZNGxYvXuxYZ7PZWLx4cZ4z9OKZPOFvvSe1oZ6S1VNygrK6iqdkdYecHnGNvKvpmjm5Wu48Vd7lzGYwy+hcpEsARETcyejRoxk0aBBt27alffv2TJw4kYyMDMco9iKu5kltqKdk9ZScoKyu4ilZzc7p9oV8WFgYVquVw4cP51l/+PBhwsPD890nPDy8WNuPGTOG0aNHO+7nXjMnUlSeOVXeeblnPtyJJx8cEZGScffdd3P06FHGjh1LcnIyLVu2ZMGCBZf0yhP35Q1/692xDS2Ip2T1lJygrK7iKVnNzOn2hfyFXef69OkDnO86N2LEiHz3iYmJYfHixYwaNcqxbuHChQV2tfP398ff39/Z0aWU8dSp8tyVpx8cEZGSMWLEiAK/D4j70996EZEr4/aFPFy+69zAgQOpXr06cXH2oyEjR46kS5cuvP322/Tq1Yu5c+eydu1aPvjgAzPfhojpPO3Mhw6OiIh4P0/5W+9JbainZPWUnKCsruIpWd0xp0cU8pfrOpeUlISPz/lx+2JjY5kzZw7PP/88zz33HPXr12f+/PmaQ15KPZ35EBERuTKe1IZ6SlZPyQnK6iqektUdc3rEPPIlTXP1ioiIu1Hb5Hz6TEVExJ0Up13yiDPyJS332IamoRMREXeR2ybp+LvzqL0XERF3Upy2XoV8Pk6ePAmgketFRMTtnDx5ktDQULNjeAW19yIi4o6K0tara30+bDYbBw8epFy5clgslqt+vtzp7Pbv36+ue06iz9T59Jm6hj5X5yutn6lhGJw8eZKIiIg848LIlXNme19afy9dSZ+pa+hzdT59pq5RGj/X4rT1OiOfDx8fH2rUqOH05w0JCSk1v4QlRZ+p8+kzdQ19rs5XGj9TnYl3Lle096Xx99LV9Jm6hj5X59Nn6hql7XMtaluvQ/oiIiIiIiIiHkSFvIiIiIiIiIgHUSFfAvz9/XnxxRfx9/c3O4rX0GfqfPpMXUOfq/PpMxV3pN9L59Nn6hr6XJ1Pn6lr6HMtnAa7ExEREREREfEgOiMvIiIiIiIi4kFUyIuIiIiIiIh4EBXyIiIiIiIiIh5EhbyIiIiIiIiIB1Eh72KTJ08mKiqKgIAAOnTowOrVq82O5NHi4uJo164d5cqVo0qVKvTp04cdO3aYHcurvPbaa1gsFkaNGmV2FI/2119/cd9991GpUiUCAwNp3rw5a9euNTuWR8vJyeGFF16gdu3aBAYGUrduXV5++WU0ZquYTW29c6mtdz219c6j9t651NYXnQp5F/riiy8YPXo0L774IuvXryc6OpoePXpw5MgRs6N5rCVLljB8+HBWrlzJwoULOXv2LDfeeCMZGRlmR/MKa9asYdq0abRo0cLsKB7t+PHjdOrUiTJlyvDTTz+xdetW3n77bSpUqGB2NI/2+uuvM2XKFCZNmsS2bdt4/fXXeeONN3j//ffNjialmNp651Nb71pq651H7b3zqa0vOk0/50IdOnSgXbt2TJo0CQCbzUZkZCSPPfYYzz77rMnpvMPRo0epUqUKS5Ys4dprrzU7jkdLT0+ndevW/Oc//+GVV16hZcuWTJw40exYHunZZ5/ljz/+4Pfffzc7ile55ZZbqFq1KtOnT3esu/POOwkMDOTTTz81MZmUZmrrXU9tvfOorXcutffOp7a+6HRG3kWysrJYt24d3bt3d6zz8fGhe/furFixwsRk3iU1NRWAihUrmpzE8w0fPpxevXrl+Z2VK/Pdd9/Rtm1b+vbtS5UqVWjVqhUffvih2bE8XmxsLIsXL2bnzp0A/PnnnyxbtoyePXuanExKK7X1JUNtvfOorXcutffOp7a+6HzNDuCtUlJSyMnJoWrVqnnWV61ale3bt5uUyrvYbDZGjRpFp06daNasmdlxPNrcuXNZv349a9asMTuKV9izZw9Tpkxh9OjRPPfcc6xZs4bHH38cPz8/Bg0aZHY8j/Xss8+SlpZGo0aNsFqt5OTk8OqrrzJgwACzo0kppbbe9dTWO4/aeudTe+98auuLToW8eKzhw4ezefNmli1bZnYUj7Z//35GjhzJwoULCQgIMDuOV7DZbLRt25YJEyYA0KpVKzZv3szUqVPVsF+FL7/8ks8++4w5c+bQtGlT4uPjGTVqFBEREfpcRbyU2nrnUFvvGmrvnU9tfdGpkHeRsLAwrFYrhw8fzrP+8OHDhIeHm5TKe4wYMYLvv/+epUuXUqNGDbPjeLR169Zx5MgRWrdu7ViXk5PD0qVLmTRpEpmZmVitVhMTep5q1arRpEmTPOsaN27Mf//7X5MSeYennnqKZ599lnvuuQeA5s2bs2/fPuLi4tS4iynU1ruW2nrnUVvvGmrvnU9tfdHpGnkX8fPzo02bNixevNixzmazsXjxYmJiYkxM5tkMw2DEiBF88803/PLLL9SuXdvsSB6vW7dubNq0ifj4eMfStm1bBgwYQHx8vBr2K9CpU6dLpkrauXMntWrVMimRdzh16hQ+PnmbLavVis1mMymRlHZq611Dbb3zqa13DbX3zqe2vuh0Rt6FRo8ezaBBg2jbti3t27dn4sSJZGRk8MADD5gdzWMNHz6cOXPm8O2331KuXDmSk5MBCA0NJTAw0OR0nqlcuXKXXHcYHBxMpUqVdD3iFXriiSeIjY1lwoQJ9OvXj9WrV/PBBx/wwQcfmB3No/Xu3ZtXX32VmjVr0rRpUzZs2MA777zDkCFDzI4mpZjaeudTW+98autdQ+2986mtLwZDXOr99983atasafj5+Rnt27c3Vq5caXYkjwbku8ycOdPsaF6lS5cuxsiRI82O4dH+7//+z2jWrJnh7+9vNGrUyPjggw/MjuTx0tLSjJEjRxo1a9Y0AgICjDp16hj/+te/jMzMTLOjSSmntt651NaXDLX1zqH23rnU1hed5pEXERERERER8SC6Rl5ERERERETEg6iQFxEREREREfEgKuRFREREREREPIgKeREREREREREPokJeRERERERExIOokBcRERERERHxICrkRURERERERDyICnkRERERERERD6JCXkRERERERMSDqJAXEY/x5JNP0qdPH7NjiIiIiIuorRcpGhXyIlKorl27MmrUKLNjABAfH0/Lli3NjiEiIuJV1NaLeB4V8iLiMf7880817iIiIl5Mbb1I0aiQFxHmzZtH8+bNCQwMpFKlSnTv3p2MjAwGDx7MkiVLeO+997BYLFgsFhITE7HZbMTFxVG7dm0CAwOJjo5m3rx5eZ6za9eujBgxghEjRhAaGkpYWBgvvPAChmFc9nXzc+DAAVJSUhyN+4kTJ+jduzedO3cmOTnZZZ+NiIiIN1BbL+JdVMiLlHKHDh2if//+DBkyhG3btvHbb79xxx13YBgG7733HjExMQwbNoxDhw5x6NAhIiMjiYuL4+OPP2bq1Kls2bKFJ554gvvuu48lS5bkee7Zs2fj6+vL6tWree+993jnnXf46KOPLvu6+YmPj6d8+fJERUWxadMm2rVrR/Xq1fn1118JDw93+eckIiLiqdTWi3gfX7MDiIi5Dh06RHZ2NnfccQe1atUCoHnz5o7H/fz8CAoKcjSgmZmZTJgwgUWLFhETEwNAnTp1WLZsGdOmTaNLly6OfSMjI3n33XexWCw0bNiQTZs28e677zq+LBT2uheLj48nOjqaOXPmMGLECF5//XWGDRvm9M9DRETE26itF/E+OiMvUspFR0fTrVs3mjdvTt++ffnwww85fvx4gdvv2rWLU6dOccMNN1C2bFnH8vHHH7N79+4823bs2BGLxeK4HxMTQ0JCAjk5OcV+3fj4eDZu3MiIESP44Ycf1LCLiIgUkdp6Ee+jQl6klLNarSxcuJCffvqJJk2a8P7779OwYUP27t2b7/bp6ekA/PDDD8THxzuWrVu3XnLtnDNfNz4+njvuuIMzZ85w4sSJYr9PERGR0kptvYj3USEvIlgsFjp16sT48ePZsGEDfn5+fPPNN4C9u11OTo5j2yZNmuDv709SUhL16tXLs0RGRuZ53lWrVuW5v3LlSurXr4/Var3s617o5MmT7Nmzh+HDhzNp0iTuuecetmzZ4uyPQURExGuprRfxLrpGXqSUW7VqFYsXL+bGG2+kSpUqrFq1iqNHj9K4cWMAoqKiWLVqFYmJiZQtW5aKFSvy5JNP8sQTT2Cz2ejcuTOpqan88ccfhISEMGjQIMdzJyUlMXr0aP7xj3+wfv163n//fd5+++0ive6F/vzzT6xWK02aNKFVq1Zs3ryZ3r17s3r1asLCwkrmgxIREfFQautFvJAhIqXa1q1bjR49ehiVK1c2/P39jQYNGhjvv/++4/EdO3YYHTt2NAIDAw3A2Lt3r2Gz2YyJEycaDRs2NMqUKWNUrlzZ6NGjh7FkyRLHfl26dDEeffRR4+GHHzZCQkKMChUqGM8995xhs9mK9LoXev/9941mzZo57mdnZxs33XSTcc011xiZmZku+mRERES8g9p6Ee9jMYwC5n8QEbkKXbt2pWXLlkycONHsKCIiIuICautFzKNr5EVEREREREQ8iAp5EREREREREQ+irvUiIiIiIiIiHkRn5EVEREREREQ8iAp5EREREREREQ+iQl5ERERERETEg6iQFxEREREREfEgKuRFREREREREPIgKeREREREREREPokJeRERERERExIOokBcRERERERHxICrkRURERERERDyICnkRERERERERD6JCXkRERERERMSD/D9Ylkyba/AsbgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xpred = model.simulate(x[0,:], C, n_steps=n-1)\n",
    "\n",
    "fig,ax = plt.subplots(1,2,figsize=(12, 3))\n",
    "\n",
    "ax[0].plot(np.linspace(0,n-1,n),x[:,0],'-bs', label='True',markersize=10)\n",
    "ax[1].plot(np.linspace(0,n-1,n),x[:,1],'-bs',markersize=10)\n",
    "ax[0].plot(np.linspace(1,n-1,n-1),xpred[:,0],'--or', label='DMDc',markersize=5)\n",
    "ax[1].plot(np.linspace(1,n-1,n-1),xpred[:,1],'--or',markersize=5)\n",
    "\n",
    "ax[0].set(ylabel=r'$x_1$', xlabel=r'steps $k$')\n",
    "ax[1].set(ylabel=r'$x_2$', xlabel=r'steps $k$')\n",
    "\n",
    "ax[0].legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
